Algorytm solwera iteracyjnego
Algorytmy solwerów iteracyjnych opisane zostały są w bardzo obszerny sposób w podręczniku Józefa Saada [1]. Prawdopodobnie pierwszym solwerem iteracyjnym była metoda rozwiązywania układu równań 4 na 4 zaproponowana przez Gaussa. Omawianie ogólnej rodziny algorytmów solwerów iteracyjnych wybiega poza ramy tego podręcznika. W module tym omówimy solwer iteracyjny mający zastosowanie w problemie projekcji w którym obszar ma dowolny kształt, nie będący regularnym kwadratem.
Niech naszym przykładowym problemem do rozwiązania będzie problem projekcji dwuwymiarowej zdefiniowany na obszarze który nie ma kształtu kwadratu (na przykład problem projekcji bitmapy na pofałdowany obszar). Wówczas nasz układ równań do rozwiązania wygląda następująco
\( \begin{bmatrix} \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{\Omega}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}u_{1,1}} \\ {\color{red}u_{2,1}} \\ \vdots \\ {\color{red}u_{N_x,N_y}}\\ \end{bmatrix} \)
\( = \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_2(y)} \\ \vdots \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} \\ \end{bmatrix} \)
gdzie \( \Omega \) nie jest już kwadratem tylko naszym ``pofałdowanym'' obszarem.
Wyprowadzenie tego układu równań przebiega w taki sam sposób jak opisano to w rozdziale pierwszym dla zadaniam projekcji bitmapy, jedyną różnicą jest fakt iż naszą projekcję wykonujemy teraz na "pofałdowany" obszar. Na przykład chcemy obliczyć i zwizualizować jak wyglądać będzie bitmapa nałożona na pofałdowany materiał (np. na zasłonę).
Oczywiście dla każdej całki po obszarze \( \Omega \) możemy dokonać zmiany zmiennych na obszar kwadratowy, i wtedy w całce znajdzie się jacobian tej transformacji
\( \int_{\Omega} {\color{blue}B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} dxdy = \int_{[0,1]^2} \color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} |Jac(x,y)| dxdy \)
Chcielibyśmy użyć do faktoryzacji naszej macierzy solwera zmiennokierunkowego, jednakże nie jest to możliwe, ponieważ nie potrafimy rozdzielić zmiennych w Jacobianie, nie wiemy jak przetworzyć \( |Jac(x,y)|=|Jac_x(x)|*|Jac_y(y)| \).
Na szczęście możemy zastosować następujący algorytm iteracyjny z preconditionerem
- Naszym celem jest rozwiązanie układu równań \( Ax=y \), gdzie \( A \) to na przykład macierz naszego problemu projekcji bitmapy na ``pofałdowany'' obszar. Rozmiar macierzy \( A \) to \( N_xN_y \) natomiast jej wyrazy dane są w naszym przykładzie przez \( {A}_{i,j}=\int_{[0,1]^2} \color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} |Jac(x,y)| dx \).
- Definiujemy nasz ``preconditioner'' czyli macierz \( M \) która powstaje z macierzy \( A \) poprzez przyjęcie jacobianu równego 1. Innymi słowy wyrazy macierzy dane są \( {M}_{i,j}=\int_{[0,1]^2}\color{blue}{B^x_{i}(x)B^y_{j}(y)} \color{red}{B^x_{k}(x)B^y_{l}(y)} dx \). Oznacza to w naszym przykładzie zastąpienie projekcji na pofałdowany obszar projekcją na płaski obszar.
- Oznaczamy \( y=Mx \)
- Wybieramy dowolny ``punkt startowy'' (dowolny wektor \( y_0 \) współczynników aproksymacji B-spline'ów o rozmiarze \( N_xN_y \)). Możemy przyjąć na przykład \( {y_0}_{i,j}= \) wartość pixela w punkcie w którym znajduje się maksimum wartości B-spline'a \( B^x_i(x) B^y_j(y) \). Ale równie dobrze mogę przyjąć \( {y_0}_{i,j}=0 \) i nasz algorytm również będzie działał poprawnie (choć pewnie wymagał będzie trochę więcej iteracji).
- Następnie obliczamy tak zwane residuum \( r= M^{-1} b - M^{-1} A M^{-1} y_0 = M^{-1} (b - A M^{-1} M y_0) \). Robimy to w następujących krokach
- Oznaczamy \( M^{-1}y_0=z_0 \) i rozwiązujemy układ równań \( Mz_0=y_0 \). Zauważmy że oznacza to rozwiązanie takiego samego układu równań jak w rozdziale pierwszym.
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}{z_0}_{1,1}} \\ {\color{red}{z_0}_{2,1}} \\ \vdots \\ {\color{red}{z_0}_{N_x,N_y}}\\ \end{bmatrix}= \) \( \begin{bmatrix} {y_0}_{1,1} \\ {y_0}_{2,1} \\ \vdots \\ {y_0}_{N_x,N_y} \end{bmatrix} \)
Możemy zrobić to za pomocą solwera zmienno-kierunkowego o liniowej złożoności obliczeniowej.
- Następnie mamy \( r=M^{-1}(b-Az_0) \) gdzie \( b-Az_0=v \) to jest wektor współczynników rozwiązania, czyli \( M^{-1}v=r \). Musimy przemnożyć macierz \( A \) (z Jacobianem) przez uzyskany wektor współczynników rozwiązania. Dostaniemy wektor
\( v=b-Az_0 = \begin{bmatrix} v_{1,1} \\ v_{2,1} \\ \cdots \\ v_{N_x,N_y} \\ \end{bmatrix} = \\ \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} - \sum_{i=1,...,N_x; j=1,...,N_y }{z_0}_{i,j} \int_{\Omega} B^x_{1,p} B^y_{1,p} B^x_{i,p}B^y_{j,p} \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y) } - \sum_{i=1,...,N_x; j=1,...,N_y} {z_0}_{i,j } \int_{\Omega} B^x_{2,p }B^y_{1,p} B^x_{i,p}B^y_{j,p} \\ \cdots \\ \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y) } - \sum_{i=1,...,N_x; j=1,...,N_y }{z_0}_{i,j} \int_{\Omega } B^x_{N_x,p}B^y_{N_y,p} B^x_{i,p}B^y_{j,p} \\ \end{bmatrix} \)
- Rozwiązujemy \( Mr=v \).
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \)
\( \begin{bmatrix} {\color{red}r_{1,1}} \\ {\color{red}r_{2,1}} \\ \vdots \\ {\color{red}r_{N_x,N_y}}\\ \end{bmatrix} \)
\( = \begin{bmatrix} \int_{\Omega} BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}{B^x_{1,p}B^y_{1,p}}{B^x_{i,p}B^y_{j,p}} \\ \int_{\Omega } BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}B^x_{2,p}B^y_{1,p}B^x_{i,p}B^y_{j,p} \\ \cdots \\ \int_{\Omega } BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} - \sum_{i=1,...,N_x;j=1,...,N_y}{z_0}_{i,j}\int_{\Omega}B^x_{N_x,p}B^y_{N_y,p}B^x_{i,p}B^y_{j,p } \\ \end{bmatrix} \)
Ponownie możemy to zrobić za pomocą solwera zmienno-kierunkowego o liniowej złożoności obliczeniowej.
Mając policzone reziduum \( r \) możemy dokonać poprawki naszego proponowanego rozwiązania
\( y_1=y_0+r \)
Warunek stopu to dostatecznie małe reziduum. Innymi słowy muszę policzyć normę z reziduum (sprawdzić jak małe jest całe reziduum, na przykład \( ||r||=\left( \sum_{i=1,...,Nx;j=1,...,Ny} |r_{i,j}|^2\right)^{0.5} \) Jeśli \( ||r||>\epsilon \) to oznaczamy \( y_0:=y_1 \) i idziemy do punktu 5.
Na końcu liczymy nasze finalne rozwiązanie (współczynniki projekcji bitmapy na pofałdowany obszar) rozwiązując układ równań \( Mx_{final}=y_{final} \). Nasze finalne rozwiązanie to \( y_{final} \).
\( \begin{bmatrix} \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{1,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{2,p}B^y_{1,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{1,p}B^y_{1,p}} & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{2,p}B^y_{1,p}} & \cdots & \int_{[0,1]^2}{B^x_{N_x,p}B^y_{N_y,p}}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} {\color{red}{x_{final}}_{1,1}} \\ {\color{red}{x_{final}}_{2,1}} \\ \vdots \\ {\color{red}{x_{final}}_{N_x,N_y}}\\ \end{bmatrix} =\\ = \begin{bmatrix} {y_{final}}_{1,1} \\ {y_{final}}_{2,1} \\ \vdots \\ {y_{final}}_{N_x,N_y} \\ \end{bmatrix} \)
Ponownie możemy to zrobić za pomocą solwera zmienno-kierunkowego o liniowej złożoności obliczeniowej.